This is a documentation of the analyses for the manuscript “(Re)Building Trust? Journals’ Open Science Badges Influence Trust in Scientists.” Currenty the data can be obtained from within this file via this link. Once the manuscript is accepted the data will be also available via the OSF.
First we had to reshape data_longest and filter the trustworthiness (METI) and relativism (TSM) items. Then we had to recode the METI items and omit the empty cases.
# A function to inverse the METI items
inv7fct <- function(x) (8-as.numeric(x))
data_psych_prop_METI <-
data_longest%>%
filter(substr(variable, 1, 3) == "tru")%>% # see codebook
mutate(item = substr(variable, 5, 9))%>%
select(session, item, treat, value)%>%
pivot_wider(names_from = item,
values_from = value)%>%
group_by(session)%>%
mutate(meas_rep = 1:n(),
treat = as.factor(treat))%>%
ungroup()%>%
mutate_at(vars(exp_1:ben_4),
list(~inv7fct(.))) # recoding 1-->5, 2-->4, ...))
data_psych_prop_tsm <-
data_longest%>%
filter(substr(variable, 1, 3) == "tsm")%>% # see codebook
mutate(item = substr(variable, 1, 5))%>%
select(session, item, treat, value)%>%
pivot_wider(names_from = item,
values_from = value)%>%
group_by(session)%>%
mutate(meas_rep = 1:n(),
treat = as.factor(treat))%>%
ungroup()%>%
mutate_at(vars(starts_with("tsm")),
list(~as.numeric(as.factor(.))))
# Join the data frames
data_psych_prop_joined <- full_join(data_psych_prop_METI, data_psych_prop_tsm)%>%
mutate(treat = as_factor(treat))
## Joining, by = c("session", "treat", "meas_rep")
# create a list of questionnaire session which are not empty
PID_list <- data_psych_prop_joined%>%
mutate(frac_NAs = rowMeans(is.na(.)))%>%
filter(frac_NAs < .8)%>%
pull(session)%>%
unique()
data_psych_prop_joined%>%
mutate(frac_NAs = rowMeans(is.na(.)))%>%
select(frac_NAs)
## # A tibble: 1,116 x 1
## frac_NAs
## <dbl>
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0.857
## 6 0.857
## 7 0.857
## 8 0.857
## 9 0
## 10 0
## # … with 1,106 more rows
# data frame with valid cases (not complete NA)
data_psych_prop <- data_psych_prop_joined%>%
filter(session %in% PID_list)
# data frames with sum scales
data_scales <- data_psych_prop%>%
mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3,
exp_4, exp_5, exp_6), na.rm = T),
Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
TSM = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), na.rm = T))
data_scales_lables <- data_psych_prop%>%
mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
Treatment = fct_recode(Treatment,
`Colored badges` = "CB",
`Control Condition` = "CC",
`Greyed out badges` = "GB"),
Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3,
exp_4, exp_5, exp_6), na.rm = T),
Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
`Topic specific multiplism` = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4),
na.rm = T))
data_scales_wide <- data_scales%>%
select(Experitse, Integrity, Benevolence,
TSM, Treatment, session)%>%
gather(Variable, Value, -session, -Treatment)%>%
mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
select(-Variable, -Treatment)%>%
spread(Variable2, Value)
data_scales_lables_wide <- data_scales_lables%>%
select(Experitse, Integrity, Benevolence,
`Topic specific multiplism`, Treatment, session)%>%
gather(Variable, Value, -session, -Treatment)%>%
mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
select(-Variable, -Treatment)%>%
spread(Variable2, Value)
# Data of treatment check
data_tch <- data_longest%>%
filter(variable %in% c("tch_1", "tch_2", "tch_3", "tch_4"))%>%
# recoding "weiß nicht[don't know]" as lowest specification of ordinal variable
mutate(value = as.ordered(value))%>%
spread(variable, value)%>%
# remove complete empty cases
filter(is.na(tch_1) == F & is.na(tch_2) == F &
is.na(tch_3) == F & is.na(tch_4) == F)%>%
group_by(session)%>%
mutate(meas_rep = 1:n())%>%
ungroup()
data_tch_n <- data_longest%>%
filter(variable %in% c("tch_1", "tch_2", "tch_3", "tch_4"))%>%
# recoding "weiß nicht[don't know]" as missing
mutate(value = as.numeric(ifelse(value == "-999", NA, value)))%>%
spread(variable, value)%>%
# remove complete empty cases
filter(is.na(tch_1) == F & is.na(tch_2) == F &
is.na(tch_3) == F & is.na(tch_4) == F)%>%
group_by(session)%>%
mutate(meas_rep = 1:n())%>%
ungroup()
length(PID_list)
## [1] 270
sum(is.na(data_psych_prop%>%
filter(meas_rep == 2)%>%
pull(exp_1)))
## [1] 13
data_longest%>%
filter(variable %in% c("age", "sex", "semester"))%>%
spread(variable, value)%>%
select(-treat)%>%
mutate(age = as.numeric(age),
semester = as.numeric(semester),
sex = as.numeric(sex))%>%
skim(.)
| Name | Piped data |
| Number of rows | 257 |
| Number of columns | 4 |
| _______________________ | |
| Column type frequency: | |
| factor | 1 |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| session | 0 | 1 | FALSE | 257 | TCT: 1, Lw_: 1, NUL: 1, 7z_: 1 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| age | 4 | 0.98 | 22.89 | 2.95 | 18 | 20 | 22 | 25 | 33 | ▇▆▅▁▁ |
| semester | 4 | 0.98 | 5.86 | 3.68 | 1 | 3 | 5 | 9 | 19 | ▇▅▅▁▁ |
| sex | 4 | 0.98 | 1.71 | 0.47 | 1 | 1 | 2 | 2 | 3 | ▃▁▇▁▁ |
data_longest%>%
filter(variable %in% c("sex"))%>%
pull(value)%>%
table(.)
## .
## 1 2 3
## 75 176 2
skewn <- function(x) DescTools::Skew(x, na.rm = T)
kurto <- function(x) DescTools::Kurt(x, na.rm = T)
maxabszvalue <- function(x) max(abs(scale(na.omit(x))))
my_skim <- skim_with(numeric = sfl(skewn, kurto, maxabszvalue))
data_scales%>%
my_skim(.)
| Name | Piped data |
| Number of rows | 540 |
| Number of columns | 26 |
| _______________________ | |
| Column type frequency: | |
| factor | 3 |
| numeric | 23 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| session | 0 | 1 | FALSE | 270 | TCT: 2, Lw_: 2, NUL: 2, 7z_: 2 |
| treat | 0 | 1 | FALSE | 3 | GB: 187, CC: 180, CB: 173 |
| Treatment | 0 | 1 | FALSE | 3 | GB: 187, CC: 180, CB: 173 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist | skewn | kurto | maxabszvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| exp_1 | 13 | 0.98 | 5.43 | 1.29 | 1 | 5.00 | 6.00 | 6.0 | 7 | ▁▁▂▃▇ | -0.87 | 0.40 | 3.42 |
| exp_2 | 13 | 0.98 | 5.44 | 1.15 | 2 | 5.00 | 6.00 | 6.0 | 7 | ▂▃▆▇▅ | -0.53 | -0.21 | 2.98 |
| exp_3 | 13 | 0.98 | 5.38 | 1.26 | 2 | 5.00 | 6.00 | 6.0 | 7 | ▂▃▆▇▅ | -0.63 | -0.15 | 2.69 |
| exp_4 | 13 | 0.98 | 5.31 | 1.32 | 1 | 4.00 | 6.00 | 6.0 | 7 | ▁▁▂▃▇ | -0.74 | 0.06 | 3.25 |
| exp_5 | 13 | 0.98 | 5.10 | 1.39 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▂▃▅▇ | -0.60 | -0.12 | 2.96 |
| exp_6 | 13 | 0.98 | 5.43 | 1.24 | 1 | 5.00 | 6.00 | 6.0 | 7 | ▁▁▂▃▇ | -0.68 | 0.02 | 3.57 |
| int_1 | 13 | 0.98 | 5.32 | 1.24 | 1 | 4.00 | 6.00 | 6.0 | 7 | ▁▁▃▃▇ | -0.49 | -0.35 | 3.48 |
| int_2 | 13 | 0.98 | 5.33 | 1.31 | 1 | 4.00 | 6.00 | 6.0 | 7 | ▁▁▃▃▇ | -0.47 | -0.50 | 3.30 |
| int_3 | 13 | 0.98 | 5.23 | 1.17 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▁▅▅▇ | -0.10 | -0.78 | 3.60 |
| int_4 | 13 | 0.98 | 5.32 | 1.22 | 1 | 4.00 | 6.00 | 6.0 | 7 | ▁▁▃▃▇ | -0.45 | -0.37 | 3.53 |
| ben_1 | 13 | 0.98 | 5.23 | 1.21 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▁▅▅▇ | -0.23 | -0.45 | 3.51 |
| ben_2 | 13 | 0.98 | 5.17 | 1.26 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▁▆▅▇ | -0.22 | -0.34 | 3.31 |
| ben_3 | 13 | 0.98 | 5.27 | 1.26 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▁▃▃▇ | -0.63 | 0.26 | 3.37 |
| ben_4 | 13 | 0.98 | 5.00 | 1.23 | 1 | 4.00 | 5.00 | 6.0 | 7 | ▁▂▆▅▇ | -0.09 | -0.56 | 3.26 |
| meas_rep | 0 | 1.00 | 1.50 | 0.50 | 1 | 1.00 | 1.50 | 2.0 | 2 | ▇▁▁▁▇ | 0.00 | -2.00 | 1.00 |
| tsm_1 | 13 | 0.98 | 1.87 | 0.87 | 1 | 1.00 | 2.00 | 2.0 | 4 | ▇▇▁▃▁ | 0.69 | -0.34 | 2.46 |
| tsm_2 | 13 | 0.98 | 2.05 | 0.95 | 1 | 1.00 | 2.00 | 3.0 | 4 | ▇▇▁▅▂ | 0.51 | -0.73 | 2.06 |
| tsm_3 | 13 | 0.98 | 2.05 | 0.97 | 1 | 1.00 | 2.00 | 3.0 | 4 | ▇▇▁▃▂ | 0.62 | -0.62 | 2.01 |
| tsm_4 | 13 | 0.98 | 2.29 | 1.01 | 1 | 1.00 | 2.00 | 3.0 | 4 | ▆▇▁▆▃ | 0.23 | -1.06 | 1.69 |
| Experitse | 13 | 0.98 | 5.35 | 1.12 | 2 | 4.67 | 5.50 | 6.0 | 7 | ▁▂▅▇▅ | -0.59 | -0.09 | 3.00 |
| Integrity | 13 | 0.98 | 5.30 | 1.08 | 2 | 4.50 | 5.50 | 6.0 | 7 | ▁▃▆▇▅ | -0.28 | -0.60 | 3.05 |
| Benevolence | 13 | 0.98 | 5.17 | 1.06 | 1 | 4.25 | 5.25 | 6.0 | 7 | ▁▁▆▇▆ | -0.16 | -0.16 | 3.93 |
| TSM | 13 | 0.98 | 2.06 | 0.67 | 1 | 1.50 | 2.00 | 2.5 | 4 | ▆▆▇▂▁ | 0.33 | -0.25 | 2.88 |
First we specified two consecutive threedimensional CFA models
# onedimensional model
cfa_meti_model_1d <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4"
cfa1d_meti_1 <- cfa(cfa_meti_model_1d, data = data_psych_prop%>%filter(meas_rep == 1))
cfa1d_meti_2 <- cfa(cfa_meti_model_1d, data = data_psych_prop%>%filter(meas_rep == 2))
cfa_meti_model_1 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int =~ int_1 + int_2 + int_3 + int_4
ben =~ ben_1 + ben_2 + ben_3 + ben_4
int_1 ~~ int_2"
cfa_meti_model_2 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int =~ int_1 + int_2 + int_3 + int_4
ben =~ ben_1 + ben_2 + ben_3 + ben_4
ben_1 ~~ ben_3
ben_1 ~~ ben_2"
cfa_meti_1 <- cfa(cfa_meti_model_1, data = data_psych_prop%>%filter(meas_rep == 1))
# define a function which prints the fit
fpf <- function(x){
fm_tmp <- fitmeasures(x)
return(cat(sprintf(
"χ^2^ = %s, _df_ = %s, CFI = %s, TLI = %s, RMSEA = %s, SRMR = %s, SRMR~between~ = %s, SRMR~within~ = %s",
round(fm_tmp[c("chisq")],3),
fm_tmp[c("df")],
round(fm_tmp[c("cfi")],3),
round(fm_tmp[c("tli")],3),
round(fm_tmp[c("rmsea")],3),
round(fm_tmp[c("srmr")],3),
round(fm_tmp[c("srmr_between")],3),
round(fm_tmp[c("srmr_within")],3))))
}
# print the fit for cfa_meti_1
fpf(cfa1d_meti_1)
χ2 = 588.932, df = 77, CFI = 0.811, TLI = 0.776, RMSEA = 0.157, SRMR = 0.084, SRMRbetween = NA, SRMRwithin = NA
fpf(cfa_meti_1)
χ2 = 194.174, df = 73, CFI = 0.955, TLI = 0.944, RMSEA = 0.078, SRMR = 0.049, SRMRbetween = NA, SRMRwithin = NA
cfa_meti_2 <- cfa(cfa_meti_model_2, data = data_psych_prop%>%filter(meas_rep == 2))
fpf(cfa1d_meti_2)
χ2 = 936.555, df = 77, CFI = 0.759, TLI = 0.715, RMSEA = 0.208, SRMR = 0.099, SRMRbetween = NA, SRMRwithin = NA
fpf(cfa_meti_2)
χ2 = 212.262, df = 72, CFI = 0.961, TLI = 0.95, RMSEA = 0.087, SRMR = 0.04, SRMRbetween = NA, SRMRwithin = NA
anova(cfa1d_meti_1, cfa_meti_1)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## cfa_meti_1 73 9738.4 9853.5 194.17
## cfa1d_meti_1 77 10125.1 10225.9 588.93 394.76 4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(cfa1d_meti_2, cfa_meti_2)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## cfa_meti_2 72 8522.8 8639.9 212.26
## cfa1d_meti_2 77 9237.1 9336.5 936.55 724.29 5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In an next step we ran a two-level CFA …
# onedimensional model
mcfa1d_meti_model <- "level: 1
meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4
level: 2
meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
int_1 + int_2 + int_3 + int_4 +
ben_1 + ben_2 + ben_3 + ben_4"
mcfa_meti_model <- "level: 1
exp_w =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int_w =~ int_1 + int_2 + int_3 + int_4
ben_w =~ ben_1 + ben_2 + ben_3 + ben_4
int_1 ~~ int_2
int_3 ~~ ben_4
ben_1 ~~ ben_2
level: 2
exp_b =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
int_b =~ int_1 + int_2 + int_3 + int_4
ben_b =~ ben_1 + ben_2 + ben_3 + ben_4"
mcfa1d_meti <- cfa(mcfa1d_meti_model, data = data_psych_prop, cluster = "session")
mcfa_meti <- cfa(mcfa_meti_model, data = data_psych_prop, cluster = "session")
fpf(mcfa1d_meti)
χ2 = 764.777, df = 154, CFI = 0.894, TLI = 0.875, RMSEA = 0.087, SRMR = 0.271, SRMRbetween = 0.146, SRMRwithin = 0.125
fpf(mcfa_meti)
χ2 = 277.759, df = 145, CFI = 0.977, TLI = 0.971, RMSEA = 0.042, SRMR = 0.172, SRMRbetween = 0.091, SRMRwithin = 0.081
anova(mcfa1d_meti, mcfa_meti)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## mcfa_meti 145 18147 18484 277.76
## mcfa1d_meti 154 18616 18915 764.78 487.02 9 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# First Measurement
semTools::reliability(cfa_meti_1)["omega",]
## exp int ben
## 0.9245381 0.8287494 0.8509624
# Second Measurement
semTools::reliability(cfa_meti_2)["omega",]
## exp int ben
## 0.9490740 0.9126052 0.8796121
First we specified two consecutive onedimensional CFA models
cfa_tsm_model <- "tsm =~ a*tsm_1 + a*tsm_2 + a*tsm_3 + a*tsm_4
tsm_2 ~~ tsm_4"
cfa_tsm_1 <- cfa(cfa_tsm_model, data = data_psych_prop%>%filter(meas_rep == 1))
fpf(cfa_tsm_1)
χ2 = 5.302, df = 4, CFI = 0.991, TLI = 0.987, RMSEA = 0.035, SRMR = 0.048, SRMRbetween = NA, SRMRwithin = NA
cfa_tsm_2 <- cfa(cfa_tsm_model, data = data_psych_prop%>%filter(meas_rep == 2))
fpf(cfa_tsm_2)
χ2 = 6.072, df = 4, CFI = 0.99, TLI = 0.985, RMSEA = 0.045, SRMR = 0.055, SRMRbetween = NA, SRMRwithin = NA
In an next step, we ran a two-level CFA …
mcfa_tsm_model <- "level: 1
tsm_w =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
tsm_2 ~~ tsm_4
level: 2
tsm_b =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
tsm_2 ~~ 0*tsm_2"
mcfa_tsm <- cfa(mcfa_tsm_model, data = data_psych_prop, cluster = "session")
fpf(mcfa_tsm)
χ2 = 4.422, df = 4, CFI = 0.999, TLI = 0.996, RMSEA = 0.014, SRMR = 0.109, SRMRbetween = 0.09, SRMRwithin = 0.019
# First Measurement
semTools::reliability(cfa_tsm_1)["omega",]
## [1] 0.5264177
# Second Measurement
semTools::reliability(cfa_tsm_2)["omega",]
## [1] 0.6459191
We specified two consecutive onedimensional CFA models
cfa_tch_model <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4
tch_1 ~~ tch_4 "
cfa_tch_1 <- cfa(cfa_tch_model,
data = data_tch_n%>%
filter(meas_rep == 1))
fpf(cfa_tch_1)
χ2 = 0.476, df = 1, CFI = 1, TLI = 1.004, RMSEA = 0, SRMR = 0.002, SRMRbetween = NA, SRMRwithin = NA
cfa_tch_2 <- cfa(cfa_tch_model, data = data_tch_n%>%filter(meas_rep == 2))
fpf(cfa_tch_2)
χ2 = 0.055, df = 1, CFI = 1, TLI = 1.018, RMSEA = 0, SRMR = 0.001, SRMRbetween = NA, SRMRwithin = NA
# First Measurement
semTools::reliability(cfa_tch_1)["alpha",]
## [1] 0.9444001
# Second Measurement
semTools::reliability(cfa_tch_2)["alpha",]
## [1] 0.9063324
tibble(`1d CFA METI 1` = fitmeasures(cfa1d_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA METI 2` = fitmeasures(cfa1d_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d CFA METI 1` = fitmeasures(cfa_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d CFA METI 2` = fitmeasures(cfa_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d MCFA METI` = fitmeasures(mcfa1d_meti)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`3d MCFA METI` = fitmeasures(mcfa_meti)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TSM 1` = fitmeasures(cfa_tsm_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TSM 2` = fitmeasures(cfa_tsm_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d MCFA TSM` = fitmeasures(mcfa_tsm)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TCH 1` = fitmeasures(cfa_tch_1)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
`1d CFA TCH 2` = fitmeasures(cfa_tch_2)[c("chisq", "df", "cfi", "tli", "rmsea",
"srmr", "srmr_between", "srmr_within", "bic", "aic")],
rownames = c("χ^2^", "_df_", "CFI", "TLI", "RMSEA", "SRMR", "SRMR~between~", "SRMR~within~", "BIC", "AIC"))%>%
column_to_rownames(var = "rownames")%>%
knitr::kable(., digits = 3)
| 1d CFA METI 1 | 1d CFA METI 2 | 3d CFA METI 1 | 3d CFA METI 2 | 1d MCFA METI | 3d MCFA METI | 1d CFA TSM 1 | 1d CFA TSM 2 | 1d MCFA TSM | 1d CFA TCH 1 | 1d CFA TCH 2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| χ2 | 588.932 | 936.555 | 194.174 | 212.262 | 764.777 | 277.759 | 5.302 | 6.072 | 4.422 | 0.476 | 0.055 |
| df | 77.000 | 77.000 | 73.000 | 72.000 | 154.000 | 145.000 | 4.000 | 4.000 | 4.000 | 1.000 | 1.000 |
| CFI | 0.811 | 0.759 | 0.955 | 0.961 | 0.894 | 0.977 | 0.991 | 0.990 | 0.999 | 1.000 | 1.000 |
| TLI | 0.776 | 0.715 | 0.944 | 0.950 | 0.875 | 0.971 | 0.987 | 0.985 | 0.996 | 1.004 | 1.018 |
| RMSEA | 0.157 | 0.208 | 0.078 | 0.087 | 0.087 | 0.042 | 0.035 | 0.045 | 0.014 | 0.000 | 0.000 |
| SRMR | 0.084 | 0.099 | 0.049 | 0.040 | 0.271 | 0.172 | 0.048 | 0.055 | 0.109 | 0.002 | 0.001 |
| SRMRbetween | NA | NA | NA | NA | 0.146 | 0.091 | NA | NA | 0.090 | NA | NA |
| SRMRwithin | NA | NA | NA | NA | 0.125 | 0.081 | NA | NA | 0.019 | NA | NA |
| BIC | 10225.876 | 9336.494 | 9853.512 | 8639.946 | 18914.759 | 18484.147 | 2791.127 | 2653.377 | 5445.587 | 2079.674 | 899.974 |
| AIC | 10125.120 | 9237.119 | 9738.362 | 8522.826 | 18616.055 | 18147.038 | 2769.537 | 2632.082 | 5360.243 | 2049.381 | 875.752 |
res_tch_data <- data_tch%>%
gather(variable, value, starts_with("tch_"))%>%
group_by(treat, variable, value)%>%
summarize(freq = n())%>%
ungroup()%>%
mutate(treat = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ treat),
value = case_when(value == "-999" ~ "don't know",
T ~ value),
variable = case_when(variable == "tch_1" ~ "item 1",
variable == "tch_2" ~ "item 2",
variable == "tch_3" ~ "item 3",
variable == "tch_4" ~ "item 4",
variable == "tch_5" ~ "item 5"),
Frequency = freq)
res_tch_plot <- ggplot(res_tch_data, aes(variable, value)) +
geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
scale_size_continuous(range = c(3,15)) +
scale_color_gradient(low = "grey95", high = "grey65") +
guides(color=guide_legend(), size = guide_legend()) +
facet_wrap(~treat) +
theme_ipsum_ps() +
ggtitle("Results of the treatment check", "Frequency per item and experimental condition") +
ylab("") +
xlab("")
res_tch_plot
# res_tch_plot_pub <- ggplot(res_tch_data, aes(variable, value)) +
# geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
# scale_size_continuous(range = c(3,15)) +
# scale_color_gradient(low = "grey95", high = "grey65") +
# guides(color=guide_legend(), size = guide_legend()) +
# facet_wrap(~treat) +
# theme_ipsum_ps() +
# ylab("") +
# xlab("")
# ggsave("res_tch_plot.svg", res_tch_plot, width = 11, height = 6)
# ggsave("Fig2.jpg", res_tch_plot_pub, width = 110, height = 50, units = "mm", dpi = 300, scale = 2.4)
res_tch_data_A <- data_tch_n%>%
filter(treat != "CC")%>%
filter(tch_1 != -999)
effsize::VD.A(tch_1 ~ treat, data = res_tch_data_A)
##
## Vargha and Delaney A
##
## A estimate: 0.8803478 (large)
data_scales_lables%>%
gather(Variable, Value, Experitse, Integrity, Benevolence)%>%
ggplot(., aes(x = Treatment, y = Value)) +
geom_violin(adjust = 1.5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
facet_wrap(~ Variable, nrow = 1) +
labs(title = "Graphical overview (Hyp 1)",
subtitle = "Violinplots and means ± 1*SD",
caption = "") +
ylim(1,7) +
hrbrthemes::theme_ipsum_ps()
## Warning: Removed 39 rows containing non-finite values (stat_ydensity).
## Warning: Removed 39 rows containing non-finite values (stat_summary).
A_GB_CC <- data_scales_lables%>%
filter(treat != "CB")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
A_CC_CB <- data_scales_lables%>%
filter(treat != "GB")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
A_GB_CB <- data_scales_lables%>%
filter(treat != "CC")%>%
mutate(treat = as.character(treat))%>%
effsize::VD.A(Integrity ~ treat, data = .)%>%
.$estimate
| GB | CC | |
|---|---|---|
| CC | A = 0.59, d = 0.32 | |
| CB | A = 0.66, d = 0.57 | A = 0.58, d = 0.29 |
plot_hyp2_1 <- ggplot(data_scales_lables,
aes(x=`Topic specific multiplism`, y = Integrity)) +
geom_jitter() +
facet_wrap(~ Treatment, nrow = 1) +
labs(title = "Graphical overview (Hyp 2/3)",
subtitle = "Jitter plot per treatment") +
hrbrthemes::theme_ipsum_ps()
plot_hyp2_1 + stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 13 rows containing non-finite values (stat_smooth).
## Warning: Removed 13 rows containing missing values (geom_point).
plot_hyp2_1 + stat_smooth(method = "lm")
## Warning: Removed 13 rows containing non-finite values (stat_smooth).
## Warning: Removed 13 rows containing missing values (geom_point).
Spearman and Kendall correlations:
r_GB <- round(cor(data_scales_lables%>%
filter(treat == "GB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "GB")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_GB <- round(cor(data_scales_lables%>%
filter(treat == "GB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "GB")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
r_CC <- round(cor(data_scales_lables%>%
filter(treat == "CC")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CC")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_CC <- round(cor(data_scales_lables%>%
filter(treat == "CC")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CC")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
r_CB <- round(cor(data_scales_lables%>%
filter(treat == "CB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CB")%>%
pull(`Topic specific multiplism`),
method = "pearson", use = "pairwise.complete.obs"), 2)
t_CB <- round(cor(data_scales_lables%>%
filter(treat == "CB")%>%
pull(Integrity),
data_scales_lables%>%
filter(treat == "CB")%>%
pull(`Topic specific multiplism`),
method = "kendall", use = "pairwise.complete.obs"), 2)
| GB | CC | CB | |
|---|---|---|---|
| \(r(integrity, multiplism)\) | -0.46 | -0.36 | -0.35 |
| \(\tau(integrity, multiplism)\) | -0.33 | -0.24 | -0.26 |
data_scales_lables%>%
mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
treat == "CB" ~ "Colored badges",
treat == "CC" ~ "Control condition",
T ~ "treat"))%>%
ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) +
geom_violin(adjust = 1.5, alpha = .5) +
stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
coord_flip() +
labs(title = "Graphical overview (Hyp 4)",
subtitle = "Violinplots and means ± 1*SD",
caption = "") +
xlab("") +
ylim(1,4) +
hrbrthemes::theme_ipsum_ps()
## Warning: Removed 13 rows containing non-finite values (stat_ydensity).
## Warning: Removed 13 rows containing non-finite values (stat_summary).
# fig4 <- data_scales_lables%>%
# mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
# treat == "CB" ~ "Colored badges",
# treat == "CC" ~ "Control condition",
# T ~ "treat"))%>%
# ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) +
# geom_violin(adjust = 1.5, alpha = .5) +
# stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
# coord_flip() +
# xlab("") +
# ylim(1,4) +
# hrbrthemes::theme_ipsum_ps()
# ggsave("Fig4.jpg", fig4, width = 120, height = 70, units = "mm", dpi = 300, scale = 1.5)
A_mult_GB_CC <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "CB")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_GB_CC <- qnorm(A_mult_GB_CC)*sqrt(2)
A_mult_CC_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "GB")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_CC_CB <- qnorm(A_mult_CC_CB)*sqrt(2)
A_mult_GB_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat,
data = data_scales_lables%>%
filter(treat != "CC")%>%
mutate(treat = as.character(treat)))$estimate
d_mult_GB_CB <- qnorm(A_mult_GB_CB)*sqrt(2)
| GB | CC | |
|---|---|---|
| CC | A = 0.43, d = -0.26 | |
| CB | A = 0.43, d = -0.25 | A = 0.5, d = 0.01 |
All of the following analyses are based on the data frame object data_scales_wide which is why we describe it here somewhat more detailed.
All analyses are based on measured variables Integrity (Integrity, is information source sincere, honest, just, unselfish and fair?) and Topic Specific Multiplism (TSM, is knowledge and knowing about a topic arbitrary?). As this data set is in wide format, the experimental conditions are encoded wtihin the variable names:
GB means, that the participants of the study could learn from the grey out badges (ans corresponding explanations) within the abstract, that the authors of the study denied to use Open PracticesCC means, that the participants of the study could not learn if or if not the authors of the study used Open PracticesCBmeans, that the participants of the study could learn from the colored badges (ans corresponding explanations) within the abstract, that the authors of the study used Open PracticesFinally, the variable session identified the study participants.
If we look descriptively at these variables:
data_scales_wide%>%
my_skim(.)
| Name | Piped data |
| Number of rows | 270 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| factor | 1 |
| numeric | 12 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| session | 0 | 1 | FALSE | 270 | TCT: 1, Lw_: 1, NUL: 1, 7z_: 1 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist | skewn | kurto | maxabszvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Benevolence_CB | 100 | 0.63 | 5.44 | 0.99 | 2.75 | 4.75 | 5.25 | 6.25 | 7 | ▁▂▇▅▆ | -0.08 | -0.71 | 2.71 |
| Benevolence_CC | 97 | 0.64 | 5.20 | 0.96 | 2.75 | 4.50 | 5.25 | 6.00 | 7 | ▁▅▇▆▃ | 0.03 | -0.61 | 2.55 |
| Benevolence_GB | 86 | 0.68 | 4.89 | 1.15 | 1.00 | 4.00 | 4.75 | 5.75 | 7 | ▁▁▇▆▅ | -0.14 | -0.02 | 3.39 |
| Experitse_CB | 100 | 0.63 | 5.63 | 1.10 | 2.33 | 4.88 | 5.83 | 6.50 | 7 | ▁▂▃▇▇ | -0.78 | 0.13 | 2.99 |
| Experitse_CC | 97 | 0.64 | 5.35 | 0.97 | 2.50 | 4.67 | 5.33 | 6.00 | 7 | ▁▃▆▇▅ | -0.25 | -0.48 | 2.94 |
| Experitse_GB | 86 | 0.68 | 5.09 | 1.19 | 2.00 | 4.29 | 5.33 | 6.00 | 7 | ▂▂▅▇▃ | -0.59 | -0.32 | 2.59 |
| Integrity_CB | 100 | 0.63 | 5.62 | 0.98 | 3.25 | 5.00 | 5.75 | 6.50 | 7 | ▂▅▇▇▇ | -0.31 | -0.83 | 2.42 |
| Integrity_CC | 97 | 0.64 | 5.34 | 0.97 | 2.75 | 4.50 | 5.50 | 6.00 | 7 | ▁▃▇▇▅ | -0.22 | -0.62 | 2.66 |
| Integrity_GB | 86 | 0.68 | 4.97 | 1.18 | 2.00 | 4.00 | 5.00 | 6.00 | 7 | ▂▅▆▇▃ | -0.09 | -0.76 | 2.53 |
| TSM_CB | 100 | 0.63 | 2.01 | 0.67 | 1.00 | 1.50 | 2.00 | 2.50 | 4 | ▇▇▇▂▁ | 0.38 | -0.34 | 2.95 |
| TSM_CC | 97 | 0.64 | 1.99 | 0.61 | 1.00 | 1.50 | 2.00 | 2.50 | 4 | ▆▆▇▁▁ | 0.25 | -0.12 | 3.30 |
| TSM_GB | 86 | 0.68 | 2.18 | 0.72 | 1.00 | 1.75 | 2.25 | 2.56 | 4 | ▅▆▇▂▂ | 0.24 | -0.46 | 2.54 |
library(naniar)
##
## Attaching package: 'naniar'
## The following object is masked from 'package:skimr':
##
## n_complete
visdat::vis_miss(data_scales_wide%>%select(-session)) +
ggtitle("Missingness per Variable") +
theme(plot.margin = margin(1, 2, 1, 1, "cm"))
Integrity_GBlibrary(patchwork)
ggplot(data_scales_wide, aes(Integrity_GB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_GB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
Integrity_CCggplot(data_scales_wide, aes(Integrity_GB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_CC)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
Integrity_CBggplot(data_scales_wide, aes(Integrity_GB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, Integrity_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, Integrity_CB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_GBggplot(data_scales_wide, aes(Integrity_GB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_GB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_GB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_CCggplot(data_scales_wide, aes(Integrity_GB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_CC)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_CC)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
TSM_CBggplot(data_scales_wide, aes(Integrity_GB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC, TSM_CB)) +
geom_miss_point(alpha = .2) +
theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB, TSM_CB)) +
geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
theme_ipsum_ps()
d_matrix_missings <- miceadds::mi_dstat(data_scales_wide%>%select(starts_with("Int"), starts_with("TSM")))%>%
round(.,4)
knitr::kable(d_matrix_missings, caption = "Boxplot of Cohen's d of missing/not-missing per variable")
| Integrity_CB | Integrity_CC | Integrity_GB | TSM_CB | TSM_CC | TSM_GB | |
|---|---|---|---|---|---|---|
| Integrity_CB | NaN | -0.2103 | -0.0096 | NaN | -0.1009 | -0.3312 |
| Integrity_CC | -0.1596 | NaN | -0.0074 | 0.1056 | NaN | 0.3622 |
| Integrity_GB | 0.1045 | 0.1696 | NaN | -0.1339 | 0.1018 | NaN |
| TSM_CB | NaN | -0.2103 | -0.0096 | NaN | -0.1009 | -0.3312 |
| TSM_CC | -0.1596 | NaN | -0.0074 | 0.1056 | NaN | 0.3622 |
| TSM_GB | 0.1045 | 0.1696 | NaN | -0.1339 | 0.1018 | NaN |
boxplot(as.vector(d_matrix_missings))
title("Boxplot of Cohen's d of missing/not-missing per variable")
M <- 1000
out <- mice(data = data_scales_wide%>%select(-session),
m = M,
meth=c("norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm"),
diagnostics = FALSE,
printFlag = F,
seed = 83851)
out_first10 <- mice(data = data_scales_wide%>%select(-session),
m = 10,
meth=c("norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm",
"norm","norm","norm"),
diagnostics = FALSE,
printFlag = F,
seed = 83851)
densityplot(out_first10)
plot(out_first10)
# Set up the matrices for the estimates ##############
# setup of matrices to store multiple estimates
mulest_hyp1 <- matrix(0,nrow=M,ncol=3)
# and covariance matrices
covwithin_hyp1 <- matrix(0,nrow=3,ncol=3)
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
within_hyp1 <- lm(cbind(Integrity_GB,Integrity_CC,Integrity_CB)~1,
data=mice::complete(out,i)) # estimate the means of the three variables
mulest_hyp1[i,]<-coef(within_hyp1)[1:3] # store these means in the matrix `mulres`
covwithin_hyp1<-covwithin_hyp1 + 1/M * vcov(within_hyp1)[1:3,1:3] # compute the averages
}
# Compute the average of the estimates ###############
estimates_hyp1 <- colMeans(mulest_hyp1)
names(estimates_hyp1) <- c("Integrity_GB","Integrity_CC","Integrity_CB")
covbetween_hyp1 <- cov(mulest_hyp1) # between covariance matrix
covariance_hyp1 <- covwithin_hyp1 + (1+1/M)*covbetween_hyp1 # total variance
# Determine the effective and real sample sizes ######
samp_hyp1 <- nrow(data_scales_wide) # real sample size
nucom_hyp1<-samp_hyp1-length(estimates_hyp1)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp1 <- (1+1/M)*(1/length(estimates_hyp1))*
sum(diag(covbetween_hyp1 %*%
MASS::ginv(covariance_hyp1))) # ... (43)
nuold_hyp1<-(M-1)/(lam_hyp1^2) # ... (44)
nuobs_hyp1<-(nucom_hyp1+1)/(nucom_hyp1+3)*nucom_hyp1*(1-lam_hyp1) # ... (46)
nu_hyp1<- nuold_hyp1*nuobs_hyp1/(nuold_hyp1+nuobs_hyp1) # ... (47)
fracmis_hyp1 <- (nu_hyp1+1)/(nu_hyp1+3)*lam_hyp1 + 2/(nu_hyp1+3) # ... (48)
neff_hyp1<-samp_hyp1-samp_hyp1*fracmis_hyp1 # = 172 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp1<-list(covariance_hyp1)
# Test the hypotheses with bain ######################
results_hyp1 <- bain(estimates_hyp1,
"Integrity_GB<Integrity_CC<Integrity_CB;
Integrity_GB=Integrity_CC=Integrity_CB;
Integrity_GB<Integrity_CC=Integrity_CB",
n = neff_hyp1, Sigma=covariance_hyp1,
group_parameters=3, joint_parameters = 0)
print(results_hyp1)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.999 0.178 5.608 4426.129 0.974 0.830
## H2 0.000 0.279 0.000 0.000 0.000 0.000
## H3 0.039 0.259 0.152 0.152 0.026 0.022
## Hu 0.148
##
## Hypotheses:
## H1: Integrity_GB<Integrity_CC<Integrity_CB
## H2: Integrity_GB=Integrity_CC=Integrity_CB
## H3: Integrity_GB<Integrity_CC=Integrity_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
summary(results_hyp1)
## Parameter n Estimate lb ub
## 1 Integrity_GB 169.8074 4.992651 4.829868 5.155435
## 2 Integrity_CC 169.8074 5.326883 5.189787 5.463978
## 3 Integrity_CB 169.8074 5.586284 5.444758 5.727809
results_hyp1$BFmatrix
## H1 H2 H3
## H1 1.000000e+00 49594808 3.700934e+01
## H2 2.016340e-08 1 7.462341e-07
## H3 2.702021e-02 1340062 1.000000e+00
path_mod <- "Integrity_GB ~ TSM_GB
Integrity_CC ~ TSM_CC
Integrity_CB ~ TSM_CB"
# Set up the matrices for the estimates ##############
best_hyp2 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp2 <- matrix(0,nrow=3,ncol=3) # and covariance matrices
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
path_fit <- sem(path_mod,
data = mice::complete(out, i),
std.lv = T) # estimate the path coefficients
best_hyp2[i,] <- parameterestimates(path_fit, standardized = T)%>% # store path coefficients
filter(op == "~")%>%
pull(std.all)
covwithin_hyp2 <- covwithin_hyp2 + # compute the average of the covariance matrices
1/M * lavInspect(path_fit,
"vcov.std.all")[c("Integrity_GB~TSM_GB",
"Integrity_CC~TSM_CC",
"Integrity_CB~TSM_CB"),
c("Integrity_GB~TSM_GB",
"Integrity_CC~TSM_CC",
"Integrity_CB~TSM_CB")]
}
# Compute the average of the estimates ###############
estimates_hyp2 <- colMeans(best_hyp2)
names(estimates_hyp2) <- c("Int_on_TSM_GB", "Int_on_TSM_CC", "Int_on_TSM_CB")
round(estimates_hyp2, 2)
## Int_on_TSM_GB Int_on_TSM_CC Int_on_TSM_CB
## -0.44 -0.41 -0.32
results_hyp2 <- bain(estimates_hyp2,
"Int_on_TSM_GB < 0 & Int_on_TSM_CC < 0 & Int_on_TSM_CB < 0;
Int_on_TSM_GB = 0 & Int_on_TSM_CC = 0 & Int_on_TSM_CB = 0",
Sigma = covariance_hyp2,
n = neff_hyp2,
group_parameters = 3,
joint_parameters = 0,
standardize = F)
print(results_hyp2)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 1.000 0.158 6.326 12953840.149 1.000 0.864
## H2 0.000 1.285 0.000 0.000 0.000 0.000
## Hu 0.136
##
## Hypotheses:
## H1: Int_on_TSM_GB<0&Int_on_TSM_CC<0&Int_on_TSM_CB<0
## H2: Int_on_TSM_GB=0&Int_on_TSM_CC=0&Int_on_TSM_CB=0
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp2$BFmatrix
## H1 H2
## H1 1.00000e+00 8.214022e+21
## H2 1.21743e-22 1.000000e+00
results_hyp3 <- bain(estimates_hyp2,
"Int_on_TSM_GB < Int_on_TSM_CC < Int_on_TSM_CB;
(Int_on_TSM_GB, Int_on_TSM_CC) < Int_on_TSM_CB;
Int_on_TSM_GB = Int_on_TSM_CC = Int_on_TSM_CB",
Sigma = covariance_hyp2,
n = neff_hyp2,
group_parameters = 3,
joint_parameters = 0,
standardize = T)
print(results_hyp3)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.541 0.177 3.058 5.488 0.122 0.118
## H2 0.820 0.354 2.319 8.332 0.093 0.089
## H3 10.086 0.513 19.645 19.645 0.785 0.755
## Hu 0.038
##
## Hypotheses:
## H1: Int_on_TSM_GB<Int_on_TSM_CC<Int_on_TSM_CB
## H2: (Int_on_TSM_GB,Int_on_TSM_CC)<Int_on_TSM_CB
## H3: Int_on_TSM_GB=Int_on_TSM_CC=Int_on_TSM_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp3$BFmatrix
## H1 H2 H3
## H1 1.0000000 1.318921 0.1556656
## H2 0.7581957 1.000000 0.1180250
## H3 6.4240278 8.472784 1.0000000
# Set up the matrices for the estimates ##############
mulest_hyp4 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp4 <- matrix(0,nrow=3,ncol=3) # and covariance matrices
# Estimate the coefficients for each data frame ######
for(i in 1:M) {
within_hyp4 <- lm(cbind( TSM_GB, TSM_CC, TSM_CB) ~ 1,
data=mice::complete(out,i)) # estimate the means of the three variables
mulest_hyp4[i,]<-coef(within_hyp4)[1:3] # store these means in the matrix `mulres`
covwithin_hyp4<-covwithin_hyp4 + 1/M * vcov(within_hyp4)[1:3,1:3] # compute the averages
}
# Compute the average of the estimates ###############
estimates_hyp4 <- colMeans(mulest_hyp4)
names(estimates_hyp4) <- c("TSM_GB","TSM_CC","TSM_CB")
covbetween_hyp4 <- cov(mulest_hyp4) # between covariance matrix
covariance_hyp4 <- covwithin_hyp4 + (1+1/M)*covbetween_hyp4 # total variance
# Determine the effective and real sample sizes ######
samp_hyp4 <- nrow(data_scales_wide) # real sample size
nucom_hyp4<-samp_hyp4-length(estimates_hyp4)
# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp4 <- (1+1/M)*(1/length(estimates_hyp4))*
sum(diag(covbetween_hyp4 %*%
MASS::ginv(covariance_hyp4))) # ... (43)
nuold_hyp4<-(M-1)/(lam_hyp4^2) # ... (44)
nuobs_hyp4<-(nucom_hyp4+1)/(nucom_hyp4+3)*nucom_hyp4*(1-lam_hyp4) # ... (46)
nu_hyp4<- nuold_hyp4*nuobs_hyp4/(nuold_hyp4+nuobs_hyp4) # ... (47)
fracmis_hyp4 <- (nu_hyp4+1)/(nu_hyp4+3)*lam_hyp4 + 2/(nu_hyp4+3) # ... (48)
neff_hyp4<-samp_hyp4-samp_hyp4*fracmis_hyp4 # = 172 approx. 2/3* 270
# coerce `covariance` to a list
covariance_hyp4<-list(covariance_hyp4)
# Test the hypotheses with bain ######################
results_hyp4 <- bain(estimates_hyp4,
"TSM_GB=TSM_CC=TSM_CB;
TSM_GB>TSM_CC>TSM_CB;
(TSM_GB,TSM_CC)>TSM_CB",
n = neff_hyp4, Sigma=covariance_hyp4,
group_parameters=3, joint_parameters = 0)
print(results_hyp4)
## Bayesian informative hypothesis testing for an object of class numeric:
##
## Fit Com BF.u BF.c PMPa PMPb
## H1 0.293 0.507 0.578 0.578 0.094 0.081
## H2 0.657 0.182 3.609 8.598 0.589 0.507
## H3 0.657 0.339 1.938 3.734 0.316 0.272
## Hu 0.140
##
## Hypotheses:
## H1: TSM_GB=TSM_CC=TSM_CB
## H2: TSM_GB>TSM_CC>TSM_CB
## H3: (TSM_GB,TSM_CC)>TSM_CB
##
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp4$BFmatrix
## H1 H2 H3
## H1 1.000000 0.1602312 0.2984187
## H2 6.240981 1.0000000 1.8624253
## H3 3.350996 0.5369343 1.0000000